MC++: A Parallel, Portable, Monte Carlo Neutron Transport Code in C++

MC++ is an implicit multi-group Monte Carlo neutron transport code written in C++ and based on the Parallel Object-Oriented Methods and Applications (POOMA) class library. MC++ runs in parallel on and is portable to a wide variety of platforms, including MPPs, SMPs, and clusters of UNIX workstations. MC++ is being developed to provide transport capabilities to the Accelerated Strategic Computing Initiative (ASCI). It is also intended to form the basis of the first transport physics framework (TPF), which is a C++ class library containing appropriate abstractions, objects, and methods for the particle transport problem. The transport problem is briefly described, as well as the current status and algorithms in MC+t for solving the transport equation. The alpha version of the POOMA class library is also discussed, along with the implementation of the transport solution algorithms using POOMA. Finally, a simple test problem is defined and performance and physics results from this problem are discussed on a variety of platforms. Description of the General Problem We begin with the time-dependent transport equation for multiplying systems (neglecting delayed neutrons)2 [Lewis & Miller 19841 1 [ l / v z + f i . 4 + X ( P , E ) a u,(P,E, t ) = Q ( P , f i , E , t ) + ' jdEjd&C;,(P, E +E, 8' fi)y(P, &, E, t ) + Where & represents the direction of travel and u,( P, fi, E, t ) represents the angular flux, such that 7 9 . (-.*-!,--? p*#.--..* . . ,., . tu..;:,., 1.. . * ; '6 . .I ~ .i .. J? IS UpJyED 1. Currently a graduate research assistant working with Stephen Lee at LANL. '1 2. Cross section notation has been modified from the reference to be self consistent in this document, which generally uses [Bell & Glasstone 19701 notation. v(P, 8, E, t)dVdQdEdt represents the total of the path lengths traveled during dt by all particles in the incremental phase space volume dVdQdE . The term Q represents external sources of neutrons, Xs( P , E + E, & 8) is the scattering cross section, and Zf( P, E ) is the fission cross section. The mean number of fission neutrons produced in a fission caused by an incident neutron of energy E is given by 2) ( E ) , and x ( E ) is the energy spectrum of fission neutrons such that v(E)Zf(P, E)x(E')dEdQ is the probable number of fission neutrons produced at P , with energies within d E about E , within the cone of angles dQ' about 8, per path length traveled by neutrons with energy E. In many transport calculations, including MC++, time-dependence is treated implicitly. Therefore, re-writing (1) in a time-independent form and assuming no external sources, A system containing fissionable material is said to be critical if there is a self-sustaining timeindependent chain reaction in the absence of external sources of neutrons. That is, after sufficient time has passed, a time-independent, asymptotic distribution of neutrons will exist such that the rate of fission neutron production is just equal to the losses due to absorption and leakage from the system. If such an equilibrium cannot be established, the neutron distribution will either increase or decrease exponentially in time. Systems are said to be subcritical if the population decreases, and supercritical if it increases. To treat the criticality problem, MC++ introduces an auxiliary eigenvalue, k. The average number of neutrons per fission event, v, is replaced by v k , and k is then adjusted to obtain a time-independent asymptotic solution to the transport equation [Bell & Glasstone 19701. Therefore, at k = keff , the effective multiplication factor, this solution is found. This amounts to varying the number of neutrons emitted per fission by I k . Therefore, we cast (2) as an eigenvalue problem and write [8 + o(P, E)]v(P, 8, E ) = j dEjd8 'Cs ( P , E + E, 81 a)v( P , 8, E ) + x ( E ) / k j d E u ( E ) Z f ( P , E) jd&v;(P , 8, E ) (3) Any chain reaction can be made critical if the number of neutrons per fission is adjusted between zero and infinity. Clearly the system is just critical if k = l . A value of kI, the number of neutrons required to make the system critical isfewer than that available in reality and hence the system is supercritical. MC++ calculates this effective multiplication factor, kef(which will be called k from now on). Description of the Algorithm All transport codes require some sort of numerical description of the geometry and composition of the system, For MCt t , the geometry description is in the form of a three-dimensional Cartesian mesh obtained from another simulation code. Each mesh element contains specific information as to its size in each dimension, location of cell vertices, material composition, and density. MC++ then stores this description (which is read from a NetCDF file) along with isotopic information supplied by the user to fully describe the materials and composition of the problem. Equation (3) is then solved in a simple iterative scheme that is widely employed in other transport codes [Briesmeister 19931 and is described briefly here. The calculation is started by guessing an initial spatial distribution of neutrons. In MC++, this initial guess is a simple scheme that places neutrons in cells containing fissile material in a “roundrobin” manner until all particles are exhausted. This, along with an initial guess for the system keffective (supplied by the user) begins the first iteration in MC++. This is called the “first generation” (or “cycle”) of our neutron population. In MCtt , k can be thought of as the ratio between the number of neutrons in successive generations, with fission events being regarded as the “birth event” that separates the generations. The mean number of fission neutrons produced in fission events are estimated and these generated neutrons are stored as the source points for the next cycle. A cycle is therefore defined as the life of all neutrons in the problem from birth (by fission) to death (by escape or capture’). Particles in the next cycle are started isotropically at the location at which the birth took place. The user controls the nominal number of particles to track per cycle and the number of cycles during which to accumulate the results. The user can also specify the number of initial cycles to “skip” to allow the neutron population to stabilize before accumulating the results. M C t t estimates the distance to collision in a given mesh cell by computing cdist = -ln(c)/Xt(E, P) (4) for each particle in the problem. Here, 5 is a uniformly distributed random number on the interval [0,1), and E, is the total macroscopic cross section in the mesh cell in which the particle resides. MC++ also computes the particle trajectory distance to the nearest boundary (which could be a mesh, problem, or material boundary) [Nolen et. al. 19961. MC++ then selects the smaller of these distances (since this is the closer event) and does the appropriate physical interaction for each particle (collision or boundary crossing event). In the Monte Carlo estimation of k, let Nbe the desired number of particles per generation (specified by the user). During the problem one tracks the particles from birth to death through a series 1. Both fission and absorption events are treated as capture events. of cycles, M. For any given fission event, the number of expected neutrons is given by in our formulation. In (5), W is the weight of the particle, and is a statistical contrivance which is a measure of an individual particle’s relative “importance” in the population. The sum of (5) over all neutron tracks in the current cycle gives the total number of particles to be tracked in the next cycle (i.e., these points are tallied and stored for the next cycle). When the next cycle begins, all particle weights are normalized by setting them equal to the ratio of N and the number of actual source points written out in the previous cycle (from (5)) . All tallies are therefore also normalized to N, rather than the sum of (5) for each cycle. Note also from (5) the effect of the v/k term in computing the number of neutrons per fission event, as shown in (3). During the cycle, MC++ accumulates information about the likelihood and result of specific events into tallies, the purpose of which is to compute estimates of k. Three different tallies, or estimators, are used in MCtt , and are the same as used in MCNP [Briesmeister 19931. The collision estimator for k is given by